﻿/* 
 * PROJECT: NyARToolkitCS
 * --------------------------------------------------------------------------------
 * This work is based on the original ARToolKit developed by
 *   Hirokazu Kato
 *   Mark Billinghurst
 *   HITLab, University of Washington, Seattle
 * http://www.hitl.washington.edu/artoolkit/
 *
 * The NyARToolkitCS is C# edition ARToolKit class library.
 * Copyright (C)2008-2009 Ryo Iizuka
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * For further information please contact.
 *	http://nyatla.jp/nyatoolkit/
 *	<airmail(at)ebony.plala.or.jp> or <nyatla(at)nyatla.jp>
 * 
 */
using System;
using System.Collections.Generic;
using System.Text;

namespace jp.nyatla.nyartoolkit.cs.core
{
    /**
     * kittlerThresholdの方式で閾値を求めます。
     * @param i_histogram
     * @return
     */
    public class NyARHistogramAnalyzer_KittlerThreshold : INyARHistogramAnalyzer_Threshold
    {
	    public int getThreshold(NyARHistogram i_histogram)
	    {
		    int i;		
		    double min=Double.MaxValue;
		    int th=0;
		    int da,sa,db,sb,pa,pb;
		    double oa,ob;
    		
		    int[] hist=i_histogram.data;
		    int n=i_histogram.length;
		    //Low側
		    da=pa=0;
		    int h;
		    for(i=0;i<n;i++){
			    h=hist[i];
			    da+=h*i;	//i*h[i]
			    pa+=h*i*i;	//i*i*h[i]
		    }
		    sa=i_histogram.total_of_data;
		    //High側(i=n-1)
		    db=0;
		    sb=0;
		    pb=0;
    		
    		
		    for(i=n-1;i>0;i--){
			    //次のヒストグラムを計算
			    int hist_count=hist[i];//h[i]
			    int hist_val =hist_count*i;  //h[i]*i
			    int hist_val2=hist_val*i;    //h[i]*i*i
			    da-=hist_val;
			    sa-=hist_count;
			    pa-=hist_val2;
			    db+=hist_val;
			    sb+=hist_count;			
			    pb+=hist_val2;

			    //初期化
			    double wa=(double)sa/(sa+sb);
			    double wb=(double)sb/(sa+sb);
			    if(wa==0 || wb==0){
				    continue;
			    }

			    oa=ob=0;
			    double ma=sa!=0?(double)da/sa:0;
			    //Σ(i-ma)^2*h[i]=Σ(i^2*h[i])+Σ(ma^2*h[i])-Σ(2*i*ma*h[i])
			    oa=((double)(pa+ma*ma*sa-2*ma*da))/sa;

			    double mb=sb!=0?(double)db/sb:0;
			    //Σ(i-mb)^2*h[i]=Σ(i^2*h[i])+Σ(mb^2*h[i])-Σ(2*i*mb*h[i])
			    ob=((double)(pb+mb*mb*sb-2*mb*db))/sb;

			    double kai=wa*Math.Log(oa/wa)+wb*Math.Log(ob/wb);
			    if(kai>0 && min>kai){
				    min=kai;
				    th=i;
			    }
			    //System.out.println(kai);

		    }
		    return th;//129//7.506713872738873
	    }

    }
}
